####################################################################################
####################################################################################
########################THE TWO FACES OF POWER-SHARING REPLICATION#############################
####################################################################################
####################################################################################

####################################################################################
########################STEP 1: PRELIMINARIES########################
####################################################################################

########1.1 DIRECTORY########
setwd("WORKING DIRECTORY")

########1.2 LIBRARIES########
library(lme4)
library(lmtest)
library(multiwayvcov)
library(stargazer)
library(plyr)
library(dplyr)
library(stargazer)
library(gplots)
library(dotwhisker)
library(broom)
library(dplyr)
library(ggrepel)

########1.3 FUNCTIONS########
get_confint<-function(model, vcovCL){
  t<-qt(.975, model$df.residual)
  ct<-coeftest(model, vcovCL)
  est<-cbind(ct[,1], ct[,1]-t*ct[,2], ct[,1]+t*ct[,2])
  colnames(est)<-c("Estimate","LowerCI","UpperCI")
  return(est)
}
cluster.se<-function(model, cluster)
{
  require(multiwayvcov)
  require(lmtest)
  vcovCL<-cluster.vcov(model, cluster)
  coef<-coeftest(model, vcovCL)
  return(coef)
}


####################################################################################
########################STEP 2: DATA########################
####################################################################################

########2.1 READ IN DATA########
twof_rep <- read.csv("two_faces_power_sharing_replication_data.csv")

########2.2 DEFINE VARIABLES########
controls_d10 <- "log(gdppc_l10) + log(pop_l10) + log(oilrents_l10 + 0.00001) + fractionalization1_l10 + nopluralitysum_l10 + postconflict_l10 + d10_victory_neg_l10 + log(d10_bdead_allc_pc_l10 + 1) + d10_un_intervention_plus_l10 + year"
controls_d5 <- "log(gdppc_l5) + log(pop_l5) + log(oilrents_l5 + 0.00001) + fractionalization1_l5 + nopluralitysum_l5 + postconflict_l5 + d10_victory_neg_l5 + log(d10_bdead_allc_pc_l5 + 1) + d10_un_intervention_plus_l5 + year"
controls_d2 <- "log(gdppc_l2) + log(pop_l2) + log(oilrents_l2 + 0.00001) + fractionalization1_l2 + nopluralitysum_l2 + postconflict_l2 + d10_victory_neg_l2 + log(d10_bdead_allc_pc_l2 + 1) + d10_un_intervention_plus_l2 + year"
controls_d1 <- "log(gdppc_l1) + log(pop_l1) + log(oilrents_l1 + 0.00001) + fractionalization1_l1 + nopluralitysum_l1 + postconflict_l1 + d10_victory_neg_l1 + log(d10_bdead_allc_pc_l1 + 1) + d10_un_intervention_plus_l1 + year"


####################################################################################
########################STEP 3: MODELS########################
####################################################################################

########3.1 CREATING DF's########
for (j in c(1,2,5,10)) {
    step1 <- paste("twof_rep_up",j," <- twof_rep",sep="")
    step2 <- paste("twof_rep_down",j," <- twof_rep",sep="")
    eval(parse(text=c(step1,step2)))
}
for (j in c(1,2,5,10)) {
  for (i in c("polityc","neldaic")) {
    step3 <- paste("twof_rep_up",j,"$",i," <- ifelse(twof_rep_up",j,"$",i," < twof_rep_up",j,"$",i,"_l",j,", twof_rep_up",j,"$",i,"_l",j,", twof_rep_up",j,"$",i,")",sep="")
    step4 <- paste("twof_rep_down",j,"$",i," <- ifelse(twof_rep_up",j,"$",i," > twof_rep_down",j,"$",i,"_l",j,", twof_rep_down",j,"$",i,"_l",j,", twof_rep_down",j,"$",i,")",sep="")
    eval(parse(text=c(step3,step4)))
  }
}

########3.2 RUNNING MODELS########
models <- c()
for (i in (c("polityc","neldaic"))) {
  for (j in c(1,2,5,10)) {
    for (k in c("", "_up", "_down")) {
      a <- gsub('eldai|olity', '',i)
      ###Main Models in different specifications
      step1 <- paste(a,"_h1_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ ps1h_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step2 <- paste(a,"_h1_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + ps1h_l",j, "+ ", get(paste("controls_d",j,sep="")), "+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step3 <- paste(a,"_h2_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ ps1h_corp_l",j,"+ ps1h_lib_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step4 <- paste(a,"_h2_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + ps1h_corp_l",j,"+ ps1h_lib_l",j, "+ ", get(paste("controls_d",j,sep="")),"+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step5 <- paste(a,"_h3_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ postconflict_l",j,"* ps1h_corp_l",j,"+ postconflict_l",j,"* ps1h_lib_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step6 <- paste(a,"_h3_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + postconflict_l",j,"* ps1h_corp_l",j,"+ postconflict_l",j,"* ps1h_lib_l",j, "+ ", get(paste("controls_d",j,sep="")),"+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      ###Robustness: no-PR measure
      step7 <- paste(a,"_h1nopr_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ ps1hnopr_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step8 <- paste(a,"_h1nopr_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + ps1hnopr_l",j, "+ ", get(paste("controls_d",j,sep="")), "+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step9 <- paste(a,"_h2nopr_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ ps1h_corpnopr_l",j,"+ ps1h_libnopr_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step10 <- paste(a,"_h2nopr_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + ps1h_corpnopr_l",j,"+ ps1h_libnopr_l",j, "+ ", get(paste("controls_d",j,sep="")),"+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step11 <- paste(a,"_h3nopr_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j," + postconflict_l",j," * ps1h_corpnopr_l",j," + postconflict_l",j,"+ ps1h_libnopr_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step12 <- paste(a,"_h3nopr_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + postconflict_l",j," * ps1h_corpnopr_l",j," + postconflict_l",j,"* ps1h_libnopr_l",j, "+ ", get(paste("controls_d",j,sep="")),"+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      ###Robustness: all components separately (to check whether different effects)
      step13 <- paste(a,"_h1comp_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ gc1_l",j, "+ pr1_l",j, "+ mv1_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step14 <- paste(a,"_h1comp_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + gc1_l",j, "+ pr1_l",j, "+ mv1_l",j, "+ ", get(paste("controls_d",j,sep="")), "+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step15 <- paste(a,"_h2comp_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j,"+ gc1_corp_l",j,"+ pr1_corp_l",j,"+ mv1_corp_l",j,"+ gc1_lib_l",j,"+ pr1_lib_l",j,"+ mv1_lib_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step16 <- paste(a,"_h2comp_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + gc1_corp_l",j,"+ pr1_corp_l",j,"+ mv1_corp_l",j,"+ gc1_lib_l",j,"+ pr1_lib_l",j,"+ mv1_lib_l",j, "+ ", get(paste("controls_d",j,sep="")),"+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step17 <- paste(a,"_h3comp_d", j, k, " <- lm(",i," - ", i, "_l",j," ~", i,"_l", j," + postconflict_l",j,"* gc1_corp_l",j," + postconflict_l",j,"* pr1_corp_l",j," + postconflict_l",j,"* mv1_corp_l",j," + postconflict_l",j,"* gc1_lib_l",j," + postconflict_l",j,"* pr1_lib_l",j," + postconflict_l",j,"* mv1_lib_l",j, "+ ", get(paste("controls_d",j,sep="")), " + ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      step18 <- paste(a,"_h3comp_fe", j, k, " <- lm(",i," ~ as.factor(gwid) + postconflict_l",j,"*  + gc1_corp_l",j," + postconflict_l",j,"*  pr1_corp_l",j," + postconflict_l",j,"*  mv1_corp_l",j," + postconflict_l",j,"*  gc1_lib_l",j," + postconflict_l",j,"*  pr1_lib_l",j," + postconflict_l",j,"*  mv1_lib_l",j, "+ ", get(paste("controls_d",j,sep="")),"+ ",i,"_meanr_l",j,",twof_rep",k,ifelse(k != '', j, ''),",na.action=na.exclude)",sep="")
      models <- c(models, paste(a,"_h1_d",j,k, sep=""))
      models <- c(models, paste(a,"_h1_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h2_d",j,k, sep=""))
      models <- c(models, paste(a,"_h2_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h3_d",j,k, sep=""))
      models <- c(models, paste(a,"_h3_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h1nopr_d",j,k, sep=""))
      models <- c(models, paste(a,"_h1nopr_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h2nopr_d",j,k, sep=""))
      models <- c(models, paste(a,"_h2nopr_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h3nopr_d",j,k, sep=""))
      models <- c(models, paste(a,"_h3nopr_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h1comp_d",j,k, sep=""))
      models <- c(models, paste(a,"_h1comp_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h2comp_d",j,k, sep=""))
      models <- c(models, paste(a,"_h2comp_fe",j,k, sep=""))
      models <- c(models, paste(a,"_h3comp_d",j,k, sep=""))
      models <- c(models, paste(a,"_h3comp_fe",j,k, sep=""))
      eval(parse(text=c(step1,step2,step3,step4,step5,step6,step7,step8,step9,step10,step11,step12,step13,step14,step15,step16,step17,step18)))
    }
  }
}

########3.3 CLUSTERING SE's########
for (i in models) {
  step1 <- paste("cse_",i,"<- data.frame(cluster.se(",i,",as.integer(twof_rep$gwid))[, 2])",sep="")
  eval(parse(text=c(step1)))
  }


####################################################################################
########################STEP 4: FIGURES AND TABLES MAIN ARTICLE########################
####################################################################################

########4.1 Table I (Overview)########
table_i <- twof_rep[c("country","year","ps1h","ps1h_corp","ps1h_lib")]
table_i$ps1h_corp_period <- as.numeric(table_i$ps1h_corp)
table_i$ps1h_lib_period <- as.numeric(table_i$ps1h_lib)
table_i_lag <- table_i
table_i_lag$year <- table_i_lag$year + 1
table_i_lag$ps1h_corp_period_lag <- table_i_lag$ps1h_corp_period
table_i_lag$ps1h_lib_period_lag <- table_i_lag$ps1h_lib_period
table_i <- left_join(table_i,table_i_lag[c("country","year","ps1h_corp_period_lag","ps1h_lib_period_lag")],by=c("country","year"))
table_i$ps1h_corp_period_lag <- ifelse(is.na(table_i$ps1h_corp_period_lag),-999,table_i$ps1h_corp_period_lag)
table_i$ps1h_lib_period_lag <- ifelse(is.na(table_i$ps1h_lib_period_lag),-999,table_i$ps1h_lib_period_lag)
table_i$ps1h_corp_period <- ifelse(is.na(table_i$ps1h_corp_period),-999,table_i$ps1h_corp_period)
table_i$ps1h_lib_period <- ifelse(is.na(table_i$ps1h_lib_period),-999,table_i$ps1h_lib_period)
table_i$pschange <- ifelse(table_i$ps1h_corp_period - table_i$ps1h_corp_period_lag > 0.1 | table_i$ps1h_corp_period - table_i$ps1h_corp_period_lag < -0.1 | table_i$ps1h_lib_period - table_i$ps1h_lib_period_lag > 0.1 | table_i$ps1h_lib_period - table_i$ps1h_lib_period_lag < -0.1,1,0)
table_i$pschange <- ifelse(is.na(table_i$pschange),1,table_i$pschange)
table_i <- table_i %>% group_by(country) %>% arrange(country,year) %>% mutate(period = cumsum(pschange))#
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(from = min(year))
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(to = max(year))
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh = min(ps1h))
table_i$min_psh <- round(as.numeric(table_i$min_psh),digits=2)
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh = max(ps1h))
table_i$max_psh <- round(as.numeric(table_i$max_psh),digits=2)
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh_corp = min(ps1h_corp))
table_i$min_psh_corp <- round(as.numeric(table_i$min_psh_corp),digits=2)
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh_corp = max(ps1h_corp))
table_i$max_psh_corp <- round(as.numeric(table_i$max_psh_corp),digits=2)
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh_lib = min(ps1h_lib))
table_i$min_psh_lib <- round(as.numeric(table_i$min_psh_lib),digits=2)
table_i <- table_i %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh_lib = max(ps1h_lib))
table_i$max_psh_lib <- round(as.numeric(table_i$max_psh_lib),digits=2)
#putting them together
table_i$time_period <- ifelse(table_i$from != table_i$to, paste(table_i$from,"-",table_i$to,sep=""),table_i$from)
table_i$horizontal_PS <- ifelse(table_i$min_psh != table_i$max_psh, paste(table_i$min_psh,"-",table_i$max_psh,sep=""),table_i$max_psh)
table_i$horizontal_PS_corp <- ifelse(table_i$min_psh_corp != table_i$max_psh_corp, paste(table_i$min_psh_corp,"-",table_i$max_psh_corp,sep=""),table_i$max_psh_corp)
table_i$horizontal_PS_lib <- ifelse(table_i$min_psh_lib != table_i$max_psh_lib, paste(table_i$min_psh_lib,"-",table_i$max_psh_lib,sep=""),table_i$max_psh_lib)
table_i <- table_i[order(-table_i$max_psh),]
table_i_all <- table_i
table_i <- unique(table_i[c("country","time_period","horizontal_PS","horizontal_PS_corp","horizontal_PS_lib")])
write.csv(table_i,file="./figures/main/tablei.csv")

########4.2 Table II (Main regressions models)########
stargazer(type="text",nc_h1_fe1,pc_h1_fe1,nc_h2_fe1,pc_h2_fe1,nc_h3_fe1,pc_h3_fe1,se=c(cse_nc_h1_fe1,cse_pc_h1_fe1,cse_pc_h2_fe1,cse_nc_h2_fe1,cse_nc_h3_fe1,cse_pc_h3_fe1), title="Table II: The effects of constitutional power-sharing on democracy (main models)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses",add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes")))

########4.3 Table III (IV First Stage)########
#Constructed in Stata

########4.3 Table IV (IV Second Stage)########
#Constructed in Stata

########4.4 Figure 1 (Coefficients)########
for (i in (c("nc_h1_fe1","pc_h1_fe1","nc_h2_fe1","pc_h2_fe1","nc_h3_fe1","pc_h3_fe1"))) {
step1 <- paste("", i, "_clustered <- cluster.se(", i, ",as.integer(twof_rep$gwid))",sep="")
step2 <- paste("", i, "_df <- data.frame(subset(tidy(", i, "_clustered, conf.int = T, conf.level = 0.9), term == 'ps1h_l1' | term == 'ps1h_corp_l1' | term == 'ps1h_lib_l1'))",sep="")
step3 <- paste("", i, "_df$conf.low <- ", i, "_df$estimate - 1.645 *", i, "_df$std.error",sep="")
step4 <- paste("", i, "_df$conf.high <- ", i, "_df$estimate + 1.645 *", i, "_df$std.error",sep="")
step5 <- paste("", i, "_df$model2 <- ", i, "_df$term",sep="")
eval(parse(text=c(step1, step2, step3, step4, step5)))
}
nc_h1_fe1_df$term <- "Power-sharing"
nc_h1_fe1_df$instrumented <- "no"
nc_h1_fe1_df$dep.var <- "Nelda"
pc_h1_fe1_df$term <- "Power-sharing"
pc_h1_fe1_df$instrumented <- "no"
pc_h1_fe1_df$dep.var <- "Polity"
nc_h2_fe1_df$term <- c("Corporate power-sharing","Liberal power-sharing")
nc_h2_fe1_df$instrumented <- "no"
nc_h2_fe1_df$dep.var <- "Nelda"
pc_h2_fe1_df$term <- c("Corporate power-sharing","Liberal power-sharing")
pc_h2_fe1_df$instrumented <- "no"
pc_h2_fe1_df$dep.var <- "Polity"
iv4.1 <- data.frame(term = c("Corporate power-sharing","Liberal power-sharing"), estimate = c(-1.721,1.393), std.error = c(0.547,0.571), instrumented = "Actors", model2 = "iv4.1", dep.var = "Nelda")
iv4.2 <- data.frame(term = c("Corporate power-sharing","Liberal power-sharing"), estimate = c(-0.675,1.651), std.error = c(0.444,0.448), instrumented = "Actors", model2 = "iv4.2", dep.var = "Polity")
iv5.1 <- data.frame(term = c("Corporate power-sharing"), estimate = c(-1.362), std.error = c(0.901), instrumented = "ZColonies", model2 = "iv5.1", dep.var = "Nelda")
iv5.2 <- data.frame(term = c("Corporate power-sharing"), estimate = c(-0.238), std.error = c(0.643), instrumented = "ZColonies", model2 = "iv5.2", dep.var = "Polity")
ivmodels <- rbind(iv4.1, iv4.2, iv5.1, iv5.2)
ivmodels$conf.low <- ivmodels$estimate - ivmodels$std.error
ivmodels$conf.high <- ivmodels$estimate + ivmodels$std.error
figure1_df <- rbind.fill(nc_h1_fe1_df, pc_h1_fe1_df, nc_h2_fe1_df, pc_h2_fe1_df, ivmodels)
figure1_df <- figure1_df[c("term","estimate","std.error","instrumented","model2","dep.var","conf.low","conf.high")]
figure1_df$model <- paste(figure1_df$instrumented, figure1_df$dep.var, sep="")
figure1_df$model <- as.factor(figure1_df$model)
figure1_df$model = factor(figure1_df$model,levels(figure1_df$model)[c(3,4,1,2,5,6)])
figure1_df$term <- as.factor(figure1_df$term)
figure1_df$term = factor(figure1_df$term,levels(figure1_df$term)[c(3,1,2)])
figure1_df$instrumented <- as.factor(figure1_df$instrumented)
figure1_df$instrumented = factor(figure1_df$instrumented,levels(figure1_df$instrumented)[c(3,1,2)])
figure1_df <- figure1_df%>% arrange((model))
figure1_df$xb <- 0
figure1_df$yb<-0
figure1 <- dwplot(figure1_df,
                       vline = geom_vline(xintercept = 0, colour = "grey60", linetype = 2),
                       dot_args = list(aes(colour = factor(dep.var))),
                       whisker_args = list(aes(colour = factor(dep.var), linetype = factor(instrumented)))) +
  theme_bw() + xlab("Coefficient estimate") + ylab("") +
  geom_line(data = figure1_df, aes(x= xb, y= yb, linetype = factor(instrumented))) +
  geom_vline(xintercept = 0, colour = "grey60", linetype = 2) +
  scale_color_manual(name = "dep. var.", values = c("black","grey60"), breaks = levels(as.factor(figure1_df$dep.var))) +
  scale_linetype_manual(name = "instrumented", values = c(3,2,1), labels = c("colonies","actors","no"), na.translate=FALSE, guide = guide_legend(reverse = TRUE)) + 
  theme(legend.position = "bottom", text=element_text(family="Times"))
ggsave(file="./figures/main/figure1.jpg", figure1, width = 17, height = 8, units ="cm", dpi = 500)


####################################################################################
########################STEP 5: FIGURES AND TABLES APPENDICES########################
####################################################################################

########5.1 Table A.I. (Descriptive Statistics)########
table_a.i <- twof_rep[c("ps1h","ps1h_corp","ps1h_lib","neldaic","neldaic_l1","polityc","polityc_l1","gdppc_l1","pop_l1","oilrents_l1","fractionalization1_l1","nopluralitysum_l1","postconflict_l1","d10_victory_neg_l1", "d10_un_intervention_plus_l1","d10_bdead_allc_pc_l1","year")]
stargazer(table_a.i,type="text",title = "Table I. Descriptive statistics")

########5.2 Tables A.V. - A.VII. (Larger time lags)########
stargazer(type="text",nc_h1_fe2,pc_h1_fe2,nc_h2_fe2,pc_h2_fe2,nc_h3_fe2,pc_h3_fe2,se=c(cse_nc_h1_fe2,cse_pc_h1_fe2,cse_pc_h2_fe2,cse_nc_h2_fe2,cse_nc_h3_fe2,cse_pc_h3_fe2), title="Table III. The effects of constitutional power-sharing on democracy (main models, lag 2)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses",add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes")))
stargazer(type="text",nc_h1_fe5,pc_h1_fe5,nc_h2_fe5,pc_h2_fe5,nc_h3_fe5,pc_h3_fe5,se=c(cse_nc_h1_fe5,cse_pc_h1_fe5,cse_pc_h2_fe5,cse_nc_h2_fe5,cse_nc_h3_fe5,cse_pc_h3_fe5), title="Table IV. The effects of constitutional power-sharing on democracy (main models, lag 5)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses",add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes")))
stargazer(type="text",nc_h1_fe10,pc_h1_fe10,nc_h2_fe10,pc_h2_fe10,nc_h3_fe10,pc_h3_fe10,se=c(cse_nc_h1_fe10,cse_pc_h1_fe10,cse_pc_h2_fe10,cse_nc_h2_fe10,cse_nc_h3_fe10,cse_pc_h3_fe10), title="Table V. The effects of constitutional power-sharing on democracy (main models, lag 10)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses",add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes")))

########5.3 Tables A.VIII. - A.XI. (First-difference specification)########
stargazer(type="text",nc_h1_d1,pc_h1_d1,nc_h2_d1,pc_h2_d1,nc_h3_d1,pc_h3_d1,se=c(cse_nc_h1_d1,cse_pc_h1_d1,cse_pc_h2_d1,cse_nc_h2_d1,cse_nc_h3_d1,cse_pc_h3_d1), title="Table VI. The effects of constitutional power-sharing on democracy (first-difference, 1-year lag)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses")
stargazer(type="text",nc_h1_d2,pc_h1_d2,nc_h2_d2,pc_h2_d2,nc_h3_d2,pc_h3_d2,se=c(cse_nc_h1_d2,cse_pc_h1_d2,cse_pc_h2_d2,cse_nc_h2_d2,cse_nc_h3_d2,cse_pc_h3_d2), title="Table VII. The effects of constitutional power-sharing on democracy (first-difference, 2-year lag)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses")
stargazer(type="text",nc_h1_d5,pc_h1_d5,nc_h2_d5,pc_h2_d5,nc_h3_d5,pc_h3_d5,se=c(cse_nc_h1_d5,cse_pc_h1_d5,cse_pc_h2_d5,cse_nc_h2_d5,cse_nc_h3_d5,cse_pc_h3_d5), title="Table VIII. The effects of constitutional power-sharing on democracy  (first-difference, 5-year lag)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses")
stargazer(type="text",nc_h1_d10,pc_h1_d10,nc_h2_d10,pc_h2_d10,nc_h3_d10,pc_h3_d10,se=c(cse_nc_h1_d10,cse_pc_h1_d10,cse_pc_h2_d10,cse_nc_h2_d10,cse_nc_h3_d10,cse_pc_h3_d10), title="Table IX. The effects of constitutional power-sharing on democracy  (first-difference, 10-year lag)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses")

########5.4 Table A.XII. (Alternate PS index without proportional representation)########
stargazer(type="text",nc_h1nopr_fe1,pc_h1nopr_fe1,nc_h2nopr_fe1,pc_h2nopr_fe1,nc_h3nopr_fe1,pc_h3nopr_fe1,se=c(cse_nc_h1nopr_fe1,cse_pc_h1nopr_fe1,cse_pc_h2nopr_fe1,cse_nc_h2nopr_fe1,cse_nc_h3nopr_fe1,cse_pc_h3nopr_fe1), title="The effects of constitutional power-sharing on democracy (alternate PS index without proportional representation)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses",add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes")))

########5.5 Table A.XIII. (disaggregation into PS components)########
stargazer(type="text",nc_h1comp_fe1,pc_h1comp_fe1,nc_h2comp_fe1,pc_h2comp_fe1,nc_h3comp_fe1,pc_h3comp_fe1,se=c(cse_nc_h1comp_fe1,cse_pc_h1comp_fe1,cse_pc_h2comp_fe1,cse_nc_h2comp_fe1,cse_nc_h3comp_fe1,cse_pc_h3comp_fe1), title="Table XI. The effects of constitutional power-sharing on democracy (disaggregation into PS components)", omit=c("factor\\(gwid"),  dep.var.labels.include=FALSE, style = "ajps", notes.append = TRUE, notes.align = "l", model.numbers = F, column.labels = c("Model 1.1 (Nelda)", "Model 1.2 (Polity)", "Model 2.1 (Nelda)", "Model 2.2 (Polity)", "Model 3.1 (Nelda)", "Model 3.2 (Polity)"), notes = "Country-clustered errors in parentheses",add.lines=list(c("Country-Fixed Effects","Yes","Yes","Yes","Yes","Yes","Yes")))

########5.6 Figure A.2 (Power-sharing and preceding democracy)########
figurea.2 <- twof_rep[c("gwid","year","ps1h","ps1h_corp","ps1h_lib","ps1h_l1","ps1h_corp_l1","ps1h_lib_l1","neldaic","neldaic_l1")]
figurea.2$change_lib <- ifelse(figurea.2$ps1h_lib > figurea.2$ps1h_lib_l1,"upwards",NA)
figurea.2$change_lib <- ifelse(figurea.2$ps1h_lib < figurea.2$ps1h_lib_l1,"downwards",figurea.2$change_lib)
figurea.2$change_lib <- ifelse(figurea.2$ps1h_lib == figurea.2$ps1h_lib_l1,"no change",figurea.2$change_lib)
figurea.2$change_lib <- as.factor(figurea.2$change_lib)
figurea.2$change_corp <- ifelse(figurea.2$ps1h_corp > figurea.2$ps1h_corp_l1,"upwards",NA)
figurea.2$change_corp <- ifelse(figurea.2$ps1h_corp < figurea.2$ps1h_corp_l1,"downwards",figurea.2$change_corp)
figurea.2$change_corp <- ifelse(figurea.2$ps1h_corp == figurea.2$ps1h_corp_l1,"no change",figurea.2$change_corp)
figurea.2$change_corp <- as.factor(figurea.2$change_corp)
plotmeans(neldaic_l1 ~ change_lib,figurea.2,connect=F,family="Times",mean.labels=F,bty="l",n.label=F,barcol="black",col="black",xlab="Countries changing their degree of\nliberal power-sharing (from t-1 to t)",ylab="Previous level of\nelectoral democracy (t-1)",cex.main=1.3, cex.lab=1.2,cex.axis=1.2)
par(family = "serif")
plotmeans(neldaic_l1 ~ change_corp,figurea.2,connect=F,family="Times",mean.labels=F,bty="l",n.label=F,barcol="black",col="black",xlab="Countries changing their degree of\ncorporate power-sharing (from t-1 to t)",ylab="Previous level of\nelectoral democracy (t-1)",cex.main=1.3, cex.lab=1.2,cex.axis=1.2)
par(family = "serif")

########5.7 Table A.XIV. Interventions and mediations########
table_a.xiv <- twof_rep[c("country","year","ps1h","ps1h_corp","ps1h_lib","d20_milit_int","d20_mediation","d20_ps1h_corp_intervener","d20_ps1h_lib_intervener","d20_ps1h_corp_mediator","d20_ps1h_lib_mediator")]
table_a.xiv_lag <- table_a.xiv
table_a.xiv_lag$year <- table_a.xiv_lag$year + 1
table_a.xiv_lag$d20_ps1h_corp_intervener_lag <- table_a.xiv_lag$d20_ps1h_corp_intervener
table_a.xiv_lag$d20_ps1h_lib_intervener_lag <- table_a.xiv_lag$d20_ps1h_lib_intervener
table_a.xiv_lag$d20_ps1h_corp_mediator_lag <- table_a.xiv_lag$d20_ps1h_corp_mediator
table_a.xiv_lag$d20_ps1h_lib_mediator_lag <- table_a.xiv_lag$d20_ps1h_lib_mediator
table_a.xiv <- left_join(table_a.xiv,table_a.xiv_lag[c("country","year","d20_ps1h_corp_intervener_lag","d20_ps1h_lib_intervener_lag","d20_ps1h_corp_mediator_lag","d20_ps1h_lib_mediator_lag")],by=c("country","year"))
table_a.xiv$d20_ps1h_corp_intervener <- ifelse(is.na(table_a.xiv$d20_ps1h_corp_intervener),-999,table_a.xiv$d20_ps1h_corp_intervener)
table_a.xiv$d20_ps1h_corp_intervener_lag <- ifelse(is.na(table_a.xiv$d20_ps1h_corp_intervener_lag),-999,table_a.xiv$d20_ps1h_corp_intervener_lag)
table_a.xiv$d20_ps1h_lib_intervener <- ifelse(is.na(table_a.xiv$d20_ps1h_lib_intervener),-999,table_a.xiv$d20_ps1h_lib_intervener)
table_a.xiv$d20_ps1h_lib_intervener_lag <- ifelse(is.na(table_a.xiv$d20_ps1h_lib_intervener_lag),-999,table_a.xiv$d20_ps1h_lib_intervener_lag)
table_a.xiv$d20_ps1h_corp_mediator <- ifelse(is.na(table_a.xiv$d20_ps1h_corp_mediator),-999,table_a.xiv$d20_ps1h_corp_mediator)
table_a.xiv$d20_ps1h_corp_mediator_lag <- ifelse(is.na(table_a.xiv$d20_ps1h_corp_mediator_lag),-999,table_a.xiv$d20_ps1h_corp_mediator_lag)
table_a.xiv$d20_ps1h_lib_mediator <- ifelse(is.na(table_a.xiv$d20_ps1h_lib_mediator),-999,table_a.xiv$d20_ps1h_lib_mediator)
table_a.xiv$d20_ps1h_lib_mediator_lag <- ifelse(is.na(table_a.xiv$d20_ps1h_lib_mediator_lag),-999,table_a.xiv$d20_ps1h_lib_mediator_lag)
table_a.xiv$ext_change <- ifelse(abs(table_a.xiv$d20_ps1h_corp_intervener - table_a.xiv$d20_ps1h_corp_intervener_lag) > 0.05 | abs(table_a.xiv$d20_ps1h_lib_intervener - table_a.xiv$d20_ps1h_lib_intervener_lag) > 0.05 | abs(table_a.xiv$d20_ps1h_corp_mediator - table_a.xiv$d20_ps1h_corp_mediator_lag) > 0.05 | abs(table_a.xiv$d20_ps1h_lib_mediator - table_a.xiv$d20_ps1h_lib_mediator_lag) > 0.05,1,0)
table_a.xiv$ext_change <- ifelse(is.na(table_a.xiv$ext_change),1,table_a.xiv$ext_change)
table_a.xiv <- table_a.xiv %>% group_by(country) %>% arrange(country,year) %>% mutate(period = cumsum(ext_change))#dito
table_a.xiv$d20_ps1h_corp_intervener <- ifelse(table_a.xiv$d20_ps1h_corp_intervener == -999, NA, table_a.xiv$d20_ps1h_corp_intervener)
table_a.xiv$d20_ps1h_corp_intervener_lag <- ifelse(table_a.xiv$d20_ps1h_corp_intervener_lag == -999, NA, table_a.xiv$d20_ps1h_corp_intervener_lag)
table_a.xiv$d20_ps1h_lib_intervener <- ifelse(table_a.xiv$d20_ps1h_lib_intervener == -999, NA, table_a.xiv$d20_ps1h_lib_intervener)
table_a.xiv$d20_ps1h_lib_intervener_lag <- ifelse(table_a.xiv$d20_ps1h_lib_intervener_lag == -999, NA, table_a.xiv$d20_ps1h_lib_intervener_lag)
table_a.xiv$d20_ps1h_corp_mediator <- ifelse(table_a.xiv$d20_ps1h_corp_mediator == -999, NA, table_a.xiv$d20_ps1h_corp_mediator)
table_a.xiv$d20_ps1h_corp_mediator_lag <- ifelse(table_a.xiv$d20_ps1h_corp_mediator_lag == -999, NA, table_a.xiv$d20_ps1h_corp_mediator_lag)
table_a.xiv$d20_ps1h_lib_mediator <- ifelse(table_a.xiv$d20_ps1h_lib_mediator == -999, NA, table_a.xiv$d20_ps1h_lib_mediator)
table_a.xiv$d20_ps1h_lib_mediator_lag <- ifelse(table_a.xiv$d20_ps1h_lib_mediator_lag == -999, NA, table_a.xiv$d20_ps1h_lib_mediator_lag)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(from = min(year))
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(to = max(year))
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh = min(ps1h))
table_a.xiv$min_psh <- round(as.numeric(table_a.xiv$min_psh),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh = max(ps1h))
table_a.xiv$max_psh <- round(as.numeric(table_a.xiv$max_psh),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh_corp = min(ps1h_corp))
table_a.xiv$min_psh_corp <- round(as.numeric(table_a.xiv$min_psh_corp),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh_corp = max(ps1h_corp))
table_a.xiv$max_psh_corp <- round(as.numeric(table_a.xiv$max_psh_corp),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_psh_lib = min(ps1h_lib))
table_a.xiv$min_psh_lib <- round(as.numeric(table_a.xiv$min_psh_lib),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_psh_lib = max(ps1h_lib))
table_a.xiv$max_psh_lib <- round(as.numeric(table_a.xiv$max_psh_lib),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_d20_ps1h_corp_intervener = min(d20_ps1h_corp_intervener,na.rm=T))
table_a.xiv$min_d20_ps1h_corp_intervener <- round(as.numeric(table_a.xiv$min_d20_ps1h_corp_intervener),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_d20_ps1h_corp_intervener = max(d20_ps1h_corp_intervener,na.rm=T))
table_a.xiv$max_d20_ps1h_corp_intervener <- round(as.numeric(table_a.xiv$max_d20_ps1h_corp_intervener),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_d20_ps1h_lib_intervener = min(d20_ps1h_lib_intervener,na.rm=T))
table_a.xiv$min_d20_ps1h_lib_intervener <- round(as.numeric(table_a.xiv$min_d20_ps1h_lib_intervener),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_d20_ps1h_lib_intervener = max(d20_ps1h_lib_intervener,na.rm=T))
table_a.xiv$max_d20_ps1h_lib_intervener <- round(as.numeric(table_a.xiv$max_d20_ps1h_lib_intervener),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_d20_ps1h_corp_mediator = min(d20_ps1h_corp_mediator,na.rm=T))
table_a.xiv$min_d20_ps1h_corp_mediator <- round(as.numeric(table_a.xiv$min_d20_ps1h_corp_mediator),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_d20_ps1h_corp_mediator = max(d20_ps1h_corp_mediator,na.rm=T))
table_a.xiv$max_d20_ps1h_corp_mediator <- round(as.numeric(table_a.xiv$max_d20_ps1h_corp_mediator),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(min_d20_ps1h_lib_mediator = min(d20_ps1h_lib_mediator,na.rm=T))
table_a.xiv$min_d20_ps1h_lib_mediator <- round(as.numeric(table_a.xiv$min_d20_ps1h_lib_mediator),digits=2)
table_a.xiv <- table_a.xiv %>% group_by(country,period) %>% arrange(country,period) %>% mutate(max_d20_ps1h_lib_mediator = max(d20_ps1h_lib_mediator,na.rm=T))
table_a.xiv$max_d20_ps1h_lib_mediator <- round(as.numeric(table_a.xiv$max_d20_ps1h_lib_mediator),digits=2)
table_a.xiv$time_period <- ifelse(table_a.xiv$from != table_a.xiv$to, paste(table_a.xiv$from,"-",table_a.xiv$to,sep=""),table_a.xiv$from)
table_a.xiv$horizontal_PS <- ifelse(table_a.xiv$min_psh != table_a.xiv$max_psh, paste(table_a.xiv$min_psh,"-",table_a.xiv$max_psh,sep=""),table_a.xiv$max_psh)
table_a.xiv$horizontal_PS_corp <- ifelse(table_a.xiv$min_psh_corp != table_a.xiv$max_psh_corp, paste(table_a.xiv$min_psh_corp,"-",table_a.xiv$max_psh_corp,sep=""),table_a.xiv$max_psh_corp)
table_a.xiv$horizontal_PS_lib <- ifelse(table_a.xiv$min_psh_lib != table_a.xiv$max_psh_lib, paste(table_a.xiv$min_psh_lib,"-",table_a.xiv$max_psh_lib,sep=""),table_a.xiv$max_psh_lib)
table_a.xiv$d20_ps1h_corp_intervener <- ifelse(table_a.xiv$min_d20_ps1h_corp_intervener != table_a.xiv$max_d20_ps1h_corp_intervener, paste(table_a.xiv$min_d20_ps1h_corp_intervener,"-",table_a.xiv$max_d20_ps1h_corp_intervener,sep=""),table_a.xiv$max_d20_ps1h_corp_intervener)
table_a.xiv$d20_ps1h_lib_intervener <- ifelse(table_a.xiv$min_d20_ps1h_lib_intervener != table_a.xiv$max_d20_ps1h_lib_intervener, paste(table_a.xiv$min_d20_ps1h_lib_intervener,"-",table_a.xiv$max_d20_ps1h_lib_intervener,sep=""),table_a.xiv$max_d20_ps1h_lib_intervener)
table_a.xiv$d20_ps1h_corp_mediator <- ifelse(table_a.xiv$min_d20_ps1h_corp_mediator != table_a.xiv$max_d20_ps1h_corp_mediator, paste(table_a.xiv$min_d20_ps1h_corp_mediator,"-",table_a.xiv$max_d20_ps1h_corp_mediator,sep=""),table_a.xiv$max_d20_ps1h_corp_mediator)
table_a.xiv$d20_ps1h_lib_mediator <- ifelse(table_a.xiv$min_d20_ps1h_lib_mediator != table_a.xiv$max_d20_ps1h_lib_mediator, paste(table_a.xiv$min_d20_ps1h_lib_mediator,"-",table_a.xiv$max_d20_ps1h_lib_mediator,sep=""),table_a.xiv$max_d20_ps1h_lib_mediator)
table_a.xiv <- table_a.xiv[order(-table_a.xiv$max_psh),]
table_a.xiv <- subset(table_a.xiv, (d20_mediation == 1 | d20_milit_int == 1) & !is.na(d20_mediation) & !is.na(d20_milit_int) & !is.na(horizontal_PS))
table_a.xiv <- unique(table_a.xiv[c("country","time_period","horizontal_PS","horizontal_PS_corp","horizontal_PS_lib","d20_ps1h_corp_intervener","d20_ps1h_lib_intervener","d20_ps1h_corp_mediator","d20_ps1h_lib_mediator","period")])
previous <- unique(table_i_all[c("country","time_period","horizontal_PS","horizontal_PS_corp","horizontal_PS_lib","period")])
previous <- previous %>% group_by(country) %>% arrange(country,period) %>% mutate(horizontal_PS_previous = na.locf(horizontal_PS,na.rm=F))
previous <- previous %>% group_by(country) %>% arrange(country,period) %>% mutate(horizontal_PS_corp_previous = na.locf(horizontal_PS_corp,na.rm=F))
previous <- previous %>% group_by(country) %>% arrange(country,period) %>% mutate(horizontal_PS_lib_previous = na.locf(horizontal_PS_lib,na.rm=F))
previous$period <- previous$period + 1
previous <- unique(previous[c("country","horizontal_PS_previous","horizontal_PS_corp_previous","horizontal_PS_lib_previous","period")])
table_a.xiv <- left_join(table_a.xiv, previous, by = c("country","period"))
table_a.xiv <- unique(table_a.xiv[c("country","time_period","d20_ps1h_corp_intervener","d20_ps1h_lib_intervener","d20_ps1h_corp_mediator","d20_ps1h_lib_mediator","horizontal_PS_corp","horizontal_PS_corp_previous","horizontal_PS_lib","horizontal_PS_lib_previous")])
write.csv(table_a.xiv,file="./figures/appendix/table_a.xiv.csv")

########5.8 Figures A.6.-A.7. Panel Matching########
library(PanelMatch)
PM.corp_nelda <- PanelMatch(lag = 4, time.id = "year", unit.id = "gwid", 
                         treatment = "ps1h_corp_binary_l1", refinement.method = "ps.match", 
                         data = twof_rep, match.missing = T, 
                         covs.formula = ~ I(ps1h_lib_l1 + log(gdppc_l1) + log(pop_l1) + log(oilrents_l1 + 0.00001) + fractionalization1_l1 + nopluralitysum_l1 + postconflict_l1 + d10_victory_neg_l1 + log(d10_bdead_allc_pc_l1 + 1) + d10_un_intervention_plus_l1 + year),
                         size.match = 10, qoi = "att" ,outcome.var = "neldaic",
                         lead = 0, forbid.treatment.reversal = FALSE)
PE.corp_nelda <- PanelEstimate(sets = PM.corp_nelda, data = twof_rep)
PM.lib_nelda <- PanelMatch(lag = 4, time.id = "year", unit.id = "gwid", 
                            treatment = "ps1h_lib_binary_l1", refinement.method = "ps.match", 
                            data = twof_rep, match.missing = T, 
                            covs.formula = ~ I(ps1h_corp_l1 + log(gdppc_l1) + log(pop_l1) + log(oilrents_l1 + 0.00001) + fractionalization1_l1 + nopluralitysum_l1 + postconflict_l1 + d10_victory_neg_l1 + log(d10_bdead_allc_pc_l1 + 1) + d10_un_intervention_plus_l1 + year),
                            size.match = 10, qoi = "att" ,outcome.var = "neldaic",
                            lead = 0, forbid.treatment.reversal = FALSE)
PE.lib_nelda <- PanelEstimate(sets = PM.lib_nelda, data = twof_rep)
PM.corp_polity <- PanelMatch(lag = 4, time.id = "year", unit.id = "gwid", 
                               treatment = "ps1h_corp_binary_l1", refinement.method = "ps.match", 
                               data = twof_rep, match.missing = T, 
                               covs.formula = ~ I(ps1h_lib_l1 + log(gdppc_l1) + log(pop_l1) + log(oilrents_l1 + 0.00001) + fractionalization1_l1 + nopluralitysum_l1 + postconflict_l1 + d10_victory_neg_l1 + log(d10_bdead_allc_pc_l1 + 1) + d10_un_intervention_plus_l1 + year),
                               size.match = 10, qoi = "att" ,outcome.var = "polityc",
                               lead = 0, forbid.treatment.reversal = FALSE)
PE.corp_polity <- PanelEstimate(sets = PM.corp_polity, 
                                  data = twof_rep)
PM.lib_polity <- PanelMatch(lag = 4, time.id = "year", unit.id = "gwid", 
                             treatment = "ps1h_lib_binary_l1", refinement.method = "ps.match", 
                             data = twof_rep, match.missing = T, 
                             covs.formula = ~ I(ps1h_corp_l1 + log(gdppc_l1) + log(pop_l1) + log(oilrents_l1 + 0.00001) + fractionalization1_l1 + nopluralitysum_l1 + postconflict_l1 + d10_victory_neg_l1 + log(d10_bdead_allc_pc_l1 + 1) + d10_un_intervention_plus_l1 + year),
                             size.match = 10, qoi = "att" ,outcome.var = "polityc",
                             lead = 0, forbid.treatment.reversal = FALSE)
PE.lib_polity <- PanelEstimate(sets = PM.lib_polity, 
                                data = twof_rep)
figure_a.6a <- DisplayTreatment(title = "Treatment distribution (dichotomized corporate power-sharing)", unit.id = "gwid",
                 time.id = "year", legend.position = "none",
                 xlab = "year",x.size = 6, ylab = "Country",hide.y.axis.label = TRUE,
                 treatment = "ps1h_corp_binary_l1", data = twof_rep)
ggsave(file="./figures/appendix/figure_a.6a.jpg", figure_a.6a, width = 17, height = 8, units ="cm", dpi = 500)
figure_a.6b <- DisplayTreatment(title = "Treatment distribution (dichotomized liberal power-sharing)", unit.id = "gwid",
                 time.id = "year", legend.position = "none",
                 xlab = "year",x.size = 6, ylab = "Country",hide.y.axis.label = TRUE,
                 treatment = "ps1h_lib_binary_l1", data = twof_rep)
ggsave(file="./figures/appendix/figure_a.6b.jpg", figure_a.6b, width = 17, height = 8, units ="cm", dpi = 500)
df1 <- data.frame(matrix(unlist(PE.corp_nelda[c(1,4)]), nrow=1, byrow=T))
df1$coefficient <- "1a: Corporate PS. (Nelda)"
df2 <- data.frame(matrix(unlist(PE.lib_nelda[c(1,4)]), nrow=1, byrow=T))
df2$coefficient <- "1b: Liberal PS. (Nelda)"
df3 <- data.frame(matrix(unlist(PE.corp_polity[c(1,4)]), nrow=1, byrow=T))
df3$coefficient <- "2a: Corporate PS. (Polity)"
df4 <- data.frame(matrix(unlist(PE.lib_polity[c(1,4)]), nrow=1, byrow=T))
df4$coefficient <- "2b: Liberal PS. (Polity)"
df <- rbind(df1,df2,df3,df4)
figure_a.7 <- ggplot(df, aes(coefficient, X1)) +
  geom_point() +
  geom_errorbar(data= df,aes(ymin= X1 - 1.645 * X2,ymax= X1 + 1.645 * X2)) + geom_hline(yintercept = 0, linetype = 2) +
  xlab("Treatment variable") + ylab("ATT") + 
  theme(text=element_text(family="Times New Roman"))
ggsave(file="./figures/appendix/figure_a.7.jpg", figure_a.7, width = 17, height = 8, units ="cm", dpi = 500)

########5.9 Figures A.8.-A.9. Case selection########
cases_plot <- cbind(twof_rep, nc_h2_fe1_fitted = fitted(nc_h2_fe1), nc_h2_fe1_resid = resid(nc_h2_fe1))
cases_plot <- cases_plot[c("gwid","country","year","systid","neldaic","nc_h2_fe1_fitted","nc_h2_fe1_resid","ps1h_corp","ps1h_lib","ps1h_corp_l1","ps1h_lib_l1")]
cases_plot <- subset(cases_plot, !is.na(nc_h2_fe1_resid) & !is.na(ps1h_corp))
cases_plot <- cases_plot %>% group_by(gwid,systid) %>% arrange(gwid,systid) %>% mutate(nc_h2_fe1_resid_systid = mean(abs(nc_h2_fe1_resid)))
cases_plot <- cases_plot %>% group_by(gwid,systid) %>% arrange(gwid,systid) %>% mutate(from = min(year))
cases_plot <- cases_plot %>% group_by(gwid,systid) %>% arrange(gwid,systid) %>% mutate(to = max(year))
cases_plot <- cases_plot %>% group_by(gwid,systid) %>% arrange(gwid,systid) %>% mutate(ps1h_corp_systid = mean(ps1h_corp))
cases_plot <- cases_plot %>% group_by(gwid,systid) %>% arrange(gwid,systid) %>% mutate(ps1h_lib_systid = mean(ps1h_lib))
cases_plot <- cases_plot %>% group_by(gwid) %>% arrange(gwid) %>% mutate(ps1h_corp_systid_change = max(abs(ps1h_corp - ps1h_corp_l1)))
cases_plot <- cases_plot %>% group_by(gwid) %>% arrange(gwid) %>% mutate(ps1h_lib_systid_change = max(abs(ps1h_lib - ps1h_lib_l1)))
cases_plot$country <- as.factor(ifelse(cases_plot$country == "Serbia / Serbia and Montenegro", "Yugoslavia", paste(cases_plot$country)))
cases_plot_corp <- subset(cases_plot,ps1h_corp_systid_change > 0)
cases_plot_lib <- subset(cases_plot,ps1h_lib_systid_change > 0)
cases_plot_corp$selected_corp <- as.factor(ifelse(cases_plot_corp$country == "Lebanon" & cases_plot_corp$systid == "798b","1","0"))
cases_plot_lib$selected_lib <- as.factor(ifelse(cases_plot_lib$country == "Nigeria" & cases_plot_lib$systid == "627","1","0"))
cases_plot_corp <- unique(cases_plot_corp[c("country","gwid","from","to","ps1h_corp_systid","nc_h2_fe1_resid_systid","ps1h_corp_systid_change","selected_corp")])
cases_plot_lib <- unique(cases_plot_lib[c("country","gwid","from","to","ps1h_lib_systid","nc_h2_fe1_resid_systid","ps1h_lib_systid_change","selected_lib")])
figure_a.8 <- ggplot(cases_plot_corp, aes(x = ps1h_corp_systid, y = nc_h2_fe1_resid_systid, colour = selected_corp)) + geom_point() +
  geom_label_repel(data = subset(cases_plot_corp, ps1h_corp_systid >=0.5), aes(family="Times", label = as.character(paste(country, " (",from,")", sep=""))),
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   size = 3,
                   segment.color = 'grey50') +
  theme_bw() + xlab("Corporate power-sharing (syst.)") + ylab("Avg. residual (syst.)") +
  theme(legend.position = "bottom", text=element_text(family="Times")) +
  scale_colour_manual(values = c("black","red"), name ="case")
figure_a.8
ggsave(file="./figures/appendix/figure_a.8.jpg", figure_a.8, width = 17, height = 6.5, units ="cm", dpi = 300)

figure_a.9 <- ggplot(cases_plot_lib, aes(x = ps1h_lib_systid, y = nc_h2_fe1_resid_systid, colour = selected_lib)) + geom_point() +
  geom_label_repel(data = subset(cases_plot_lib, ps1h_lib_systid >=0.45), aes(family="Times", label = as.character(paste(country, " (",from,")", sep=""))),
                   box.padding   = 0.35, 
                   point.padding = 0.5,
                   size = 3,
                   segment.color = 'grey50') +
  theme_bw() + xlab("liberal power-sharing (syst.)") + ylab("Avg. residual (syst.)") +
  theme(legend.position = "bottom", text=element_text(family="Times")) +
  scale_colour_manual(values = c("black","red"), name ="case")
figure_a.9
ggsave(file="./figures/appendix/figure_a.9.jpg", figure_a.9, width = 17, height = 6.5, units ="cm", dpi = 300)
